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Summary 


This paper is concerned with two important elements in the high-order accurate spatial discretization 
of finite volume equations over arbitrary grids. One element is the integration of basis functions over 
arbitrary domains, which is used in expressing various spatial integrals in terms of discrete unknowns. The 
other consists of quadrature approximations to those integrals. Only polynomial basis functions applied 
to polyhedral and polygonal grids are treated here. Non- triangular polygonal faces are subdivided into a 
union of planar triangular facets, and the resulting triangulated polyhedron is subdivided into a union of 
tetrahedra. The straight line segment, triangle, and tetrahedron are thus the fundamental shapes that are 
the building blocks for all integrations and quadrature approximations. Integrals of products up to the fifth 
order are derived in a unified manner for the three fundamental shapes in terms of the position vectors of 
vertices. Results are given both in terms of tensor products and products of Cartesian coordinates. The exact 
polynomial integrals are used to obtain symmetric quadrature approximations of any degree of precision up to 
five for arbitrary integrals over the three fundamental domains. Using a coordinate-free formulation, simple 
and rational procedures are developed to derive virtually all quadrature formulas, including some previously 
unpublished. Four symmetry groups of quadrature points are introduced to derive Gauss formulas, while 
their limiting forms are used to derive Lobatto formulas. Representative Gauss and Lobatto formulas are 
tabulated. The relative efficiency of their application to polyhedral and polygonal grids is detailed. The 
extension to higher degrees of precision is discussed. 


I. Introduction 


Two of the current themes in computational physics are high-order accurate methods and the employ- 
ment of arbitrary grids. We are interested in using these ideas in the discretization of integral conservation 
laws. As applied to a finite domain, these state that the rate of change of a conserved quantity inside a 
domain is equal to integrated effects over its boundary, and possibly its net creation inside the domain. For 
three-dimensional problems, the domain is a volume element, and the boundary is a closed surface. For 
two-dimensional problems, the dimensions of the geometric elements are reduced by one. 


A general procedure for obtaining a high-order accurate discretization of equations on a given arbitrary 
grid would involve the following steps: 

1. A set of basis functions is chosen to represent the solution in some local region. This region is usually 
the finite conservation domain, but it could be centered with respect to part of the boundary of the con- 
servation domain. The basis functions are normally chosen as polynomials, but other readily integrable 
or differentiable functions can also be used. 

2. A set of discrete unknowms in the neighborhood of the representation region is chosen to reconstruct 
the local representation of the solution. These can be symmetrically located or directionally biased. 
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3. The coefficients in the representation are evaluated in terms of the conservative unknowns. This requires 
the integration of the basis functions over the conservation domain. In the most general case, a least- 
squares algebraic problem with constraints must be solved. 

4. The boundary terms, and possible source terms, are evaluated for the conservation domain. If the bound- 
ary functions are linear, and the boundary is analytically defined, the evaluation could be performed 
in closed form, provided that the basis functions are integrable. In general, a high-order quadrature in 
terms of point values is required. The solution at a boundary point could involve spatial derivatives 
(requiring the differentiation of the basis functions) if transport terms are present. It will also generally 
involve non-linear functions of the states on the two sides of the boundary. 

The present paper is restricted to a discussion of two important elements in such a procedure. They are the 
integration of the basis functions over arbitrary domains and quadrature approximations to general integrals 
over such domains. 

An arbitrary grid is defined by specifying a set of points and the lines connecting those points. For three 
dimensions, the precise surfaces that form cell boundaries must also be specified. These lines and surfaces 
should be defined so as to make the integration of basis functions and quadrature approximations easily 
performed. (An exception could be for lines and surfaces at specified global domain boundaries, where more 
complex evaluations for certain analytical shapes may be tolerated.) One therefore normally connects the 
points by straight Unes, although piecewise straight lines may be necessary when employing dual grids. In 
two dimensions, this produces plane polygonal domains in general. For three dimensions, the finite domains 
are arbitrary polyhedra, with polygonal faces of different type. For other than triangles, the faces are in 
general non-planar. If the surface of a quadrilateral face is defined as a ruled surface, the integrations of 
certain basis functions can be performed analytically, but the resulting expressions are extremely complex. 
Also, quadrature approximations would be very difficult to obtain. It is therefore best to subdivide each 
non-triangular face into a union of planar triangular facets. Efficient volume integrations can be obtained 
by subdividing the resulting triangulated polyhedron into a union of tetrahedra. The face triangulations 
define the shapes of the boundary surfaces, and therefore affect the answers for the surface and volume 
integrations. However, once the shapes of the boundary surfaces are defined, the subdivision into tetrahedra 
is just a matter of convenience in performing the volume integrations, and would not affect the final answers. 
The most efficient subdivision of a polygonal face, in terms of minimizing the number of triangular facets, 
can always be accomplished without introducing any additional interior points. But an efficient subdivision 
of a triangulated polyhedron into tetrahedra may require the introduction of interior vertices or interior 
edges. The straight line segment, triangle, and tetrahedron are thus the fundamental shapes that are the 
building blocks for all integrations and quadrature approximations. 

The simplest basis functions are polynomials, and this paper will be restricted to their integration over 
the fundamental domains. Other basis functions can be more appropriate for certain classes of problems, and 
will be considered in a future publication. A coordinate-free formulation using vector and tensor notation is 
employed, thus simplifying the derivations, and permitting a unified presentation for the three fundamental 
shapes. In this formulation, any function is expressed as a generalized Taylor series in terms of tensor products 
of the position vector, whose origin can be arbitrarily defined. The integration of a general polynomial thus 
involves the integration of tensor products of various order over the fundamental domains. It is possible to 
express the answers for all three fundamental shapes by a unified formula for a tensor product of a given 
order. Such expressions axe given for products up to the fifth order, and higher-order formulas can be easily 
derived by the same procedure. When the order of the tensor product is greater than the number of vertices 
defining the fundamental shape, the expression is no longer unique. Some terms in the unified formulas 
can then be written in terms of others. We thus obtain simplified formulas for the line segment starting 
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with the third order, the triangle for the fourth order, and the tetrahedron for the fifth order. For practical 
applications, expressions for generic monomials in terms of Cartesian coordinates are also given. From these 
one can easily obtain the expression for any specific monomial. 

Symmetric quadrature approximations to integrals over an n-dimensional simplex, including triangles 
and tetrahedra, have received much attention in recent years. A compilation of formulas is found in the 
book by Stroud (ref. 1), with references to the original papers given therein. Some additional formulas 
are found in Grundmann and Moller (ref. 2). Formulas for triangles were also given by Cowper (ref. 3) 
and Lyness and Jespersen (ref. 4), and for tetrahedra by Yu (ref. 5) and Keast (ref. 6). Some of the 
derivations of these formulas involved algebraic methods, based on roots of orthogonal polynomials, and 
were often of an ad hoc nature. Using our coordinate- free formulations, we develop very simple and rational 
procedures to derive all of the above formulas, and also some useful ones that are new to our knowledge. 
The exact polynomial integrals are used to obtain quadrature approximations of various degrees of precision 
for arbitrary integrals over the fundamental domains. For isolated fundamental shapes, the most efficient 
formulas, minimizing the number of quadrature points, are those in which the locations of the quadrature 
points are unspecified. By analogy with the one-dimensional cases, the formulas are referred as Gauss, even 
for triangles and tetrahedra. For the latter two shapes, due to the non-linearity of the equations, there 
may be more than one Gauss formula for a given degree of precision, and some could possess properties not 
shared by the one-dimensional formulas. Quadrature points could be located outside or on the boundaries, 
and some of the coefficients or weights could be negative, which may be unsuitable in some applications. 
Under some conditions, it may be desirable to employ formulas that involve one or more free parameters. 
While they are less efficient for isolated fundamental shapes, the parameters can be chosen so that all the 
weights are positive, if this is necessary. A more import ant role for the free parameter (s) is to place some 
or all the quadrature points on shape boundaries (vertices, edges, or faces). Such formulas are referred as 
Lobatto, again by analogy with the one-dimensional case. When the original finite domains are part of a grid 
consisting of arbitrary polyhedra or polygons, these boundary elements are shared by several fundamental 
shapes. Depending on the type of the polyhedra or polygons, the degree of precision, the amount of storage 
available, and the parallelization of a code, a Lobatto formula may be more efficient than a Gauss formula. 
We therefore present both Gauss and Lobatto formulas of various degrees of precision. Note that besides 
the evaluations of boundary terms and possible source terms, quadrature approximations are useful when 
calculating integrals over the conservation domain for those basis functions that are not readily integrable. 
The type of term being evaluated w r ill also play a role in choosing between a Gauss or Lobatto formula. 


II. Exact Integrals of Polynomials for the Fundamental Shapes 


Let r denote the position vector of a point in space with respect to an arbitrary origin. Using tensor 
notation, an arbitrary function /( r) can be expanded about the origin in terms of tensor products as a 
generalized Taylor series 

/( r) = /( 0) + r • (V/)o + | rr : (VV/) 0 + | rrr : (VVV/)o + ■ • • • ( 1 ) 

In our procedure for obtaining a high-order accurate discretization of equations, steps (3) and (4) require 
the integration of /( r) over the conservation domain and its boundaries, in which the evaluations of the 
integrals of the tensor products of r of various order over the fundamental domains are required. Let n be 
the number of vertices defining the fundamental shape. The line segment (n = 2) is defined by the points 
ri and r 2 > and has the length 

L = |r 2 — ri |. ( 2a ) 
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Similarly, the triangle (n = 3) is defined by the points ri, r 2 , and r 3 , and has the area 

5= ||(r 2 -ri) x (r 3 - r 2 )|, 

and the tetrahedron (n — 4) is defined by the points ri, r 2 , r 3 , and r 4 , and has the volume 


y = gl( r 2 - ri) X (r 3 - Ti) • (r 4 - ri)|. 


(26) 


(2c) 


In order to derive the formulas for the integrals of the tensor products, instead of employing the conven- 
tional simplex coordinates, it is simpler for our purposes to use the independent parametric representations 
for the three fundamental shapes shown in figure l(a-c). Here the parameters u , t;, and w range from 0 to 
1 . A point on the line segment is given by 


r = ri(l - u) + V 2 U, 


and the differential element dL is 


dL = 


dr 

du 


du = L du . 


A point on the triangular face is given by 


(за) 

( зб ) 


r = n (1 - u) + r 2 tt(l - v) + r 3 utf, 


(4a) 


and the differential element is 


dS 


dr dr 
du X dv 


du dv = 2Su du dv. 


(46) 


Similarly, a point in the tetrahedron is given by 


r = ri(l — u) + r 2 t/(l — v) -f r 3 uu(l — w) + r ^uvw, 


(5a) 


and the differential element is 


dV = 


dr dr 
du dv 


dr 

dw 


du dv dw = 6V u 2 v du dv dw. 


(56) 


Using the above changes of variables, the spatial integrals can be transformed to integrals in the parameter 
space. The derivations can be simplified considerably by making use of certain symmetries. For example, in 
evaluating / rrr dV, since the final result must be symmetric in the four vertices, the coefficients of ririri 
and T 2 r 2 r 2 must be the same, etc. Similarly, since products of powers of u, v and w commute, it follows 
that the coefficients of ririr 2 and ri^ri must be the same, etc. As a result, the integrations only involve 
a few independent integrals in the parameter space. The formulas for the exact integrals of tensor products 
over the three fundamental domains can be expressed in a unified manner in terms of tensor products of the 
position vectors of the vertices. The final results up to the fifth order are 


/ 

/ 


dL 
1 dS 
dV 

dL 
r dS 
dV 


( L\ ( |r 2 — ri| \ 

5 = i|(r 2 -r x ) x (rs-ri)l 

\V ) \f|( r 2 -ri) x (r 3 -ri)-(r 4 — ri)|/ 



(ба) 

( бб ) 
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dL 

rrr dS = 
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dV ■ -/ i j k i j 

dL j 

rrrr dS 1 — , , _ x < / / 

n(n + l)(n -f 2)(n 4- 3) i j k i 


C EE r i r ;j rfc + ^^(r i r i r j + ---)+2E] r t r « r *] 5 


(6c) 

(6d) 


EEEE r iri r fc r( + E E E( r * r * r ' rfc + ’ ' ') 

i j k 


+ E ^( T i T i T J r i + •••)+ 2 53 Efo r<r ‘ r j +■■•) + 6 E r iri r t r 2 ] ( S J (6e) 

i j i 3 * \ ^ ) 


/ 


dL j 

rrrrr d 5 n ( n + 1 )( n + 2 )( n + 3)( n + 4)iZ-^Z- fc ; m 


EEEEE rir J r * r ' r " 


+ E E E J2( Tir ' r j TkTi + ■ ■ •) + 23 E E( r * rir ^ rfc + • • ■) 

i j k l i j k 

+ 2 e E Efcw r * +•••)+ 2 E E( r < r * r w + • ■ ■) 


i j k 


* J 


+ 6 ED r^r^r, 4 )+24Er i r<r i r i r t ]( S 


( 6 /) 


Z J 


Here the range for each summation is from 1 to n. In eqs. (6d) to (6f), we have only indicated one term 
in some of the summations, the others being obtained from the independent permutations of the indices. 
For example, there are three terms in the second summation in eq. (6d), namely, r^r,, r^r*, and r^r*^. 
Similarly, there are six, three, and four terms in the second to fourth summations in eq. (6e), and ten, fifteen, 
ten, ten, and five terms in the second to sixth summations in eq. (6f), respectively. The above formulas 
can be simplified when the order of the tensor products is greater than n. It can then be shown that the 
contribution of all the summations involving an even number of indices is equal to that of the summations 
involving an odd number of indices. Thus for the line segment (n = 2), eq. (6d) can be replaced by 

n(n + lXn+2 ) E E E ^ + 2 E 


(7a) 


i j k i 

For the line segment ( n = 2) and the triangle (n = 3), eq. (6e) can be replaced by 
dL 2 


/ 


rrrr 


dS n(n + l)(n -f 2 )(n -f 3) 


i j k 

Finally, for all three shapes, eq. (6f) can be replaced by 

dL 


E E E( r *w* + • • •) + 6E r » r « r * r *] ( 


(7b) 


/ 


rrrrr dS = 


EEEEE r ^ rtrirm 


jy n(n + l)(n + 2)(n 4* 3)(n 4- 4) *-r* " fc x 

+EID'w jT k 4- • • •) 4- 2 E E E( r * r ‘ r ^ rfc + •••)+ 24 E Wr^r*] ( 5 


(7c) 


i j k 


i j k 


In practical calculations, the position vector is usually expressed in terms of Cartesian coordinates. Let 

< >=E( )< (8) 



denote the sum over the vertices of the fundamental shape. For example, < xy >= X\y\ + £2^2 for a line 
segment. From eqs. (6) and (7), we can easily derive the following relations for the integrals of generic 
monomials in Cartesian coordinates: 


/ 


dL 


1 


x dS = — < x > 
dV n 

dL 


xy dS = 


/ Ulj 

xyz dS = 


dV n(n + 1) 
dL 


S 

v t 

[< x >< y > + < xy >] 

[< x >< y >< z > 


L 

S 

V 


(9a) 

(96) 


dV n(n + l)(n + 2) 

+ < i >< yz > + < y X xz > + < z X xy > +2 < xyz >] 


' L ' 
S 
V , 


(9c) 


I 


dL 

x 2 yz dS = 


d y n(n + l)(n + 2)(n + 3) 


[(< x > 2 + < x 2 >)(< y >< z > + < yz >) 


+ 2 < y > (< x >< xz > + < x 2 z >) + 2 < z > (< x >< xy > + < x 2 y >) 

( L > 

+ 2 < xy >< xz > +4 < x >< xyz > +6 < x 2 yz >] I S j . (9 d) 

w 


When the order of the monomial is greater than n, the relations can be simplified, and one obtains 
2 


/ 


xyz dL — 


n(n 4- l)(n + 2) 


[< x >< y >< z > +2 < xyz >] L 


/ 

I 


2 dL 
x yz , c = 


[< x > 2 < yz > + < x 2 >< y >< z > 


dS n(n + l)(n + 2)(n + 3) 

+ 2 < x > (< y >< xz > + < z >< xy >) +6 < x 2 yz >] ^5^ 


(10a) 


(106) 


dL 

x 3 yz dS = 


[< x > < y >< z > 


d y n(n + l)(n + 2)(n + 3)(n + 4) 

+ 6 < x > (< xy X xz > + < x >< xyz > + < y >< x 2 z > + < z >< x 2 y >) 

+ 3(< x >< x 2 >< yz > + < y >< x 2 >< xz > + < z >< x 2 X xy >) 

L ' 

+ 12 < x 3 X y X z > +24 < x 3 yz >] ( S 


(10c) 


/■ 


dL 


x 2 y 2 z dS 


:< x > 2 < y > 2 < 2 > + < x 2 X y 2 X z > 


d y n(n + l)(n + 2)(n + 3)(n + 4) 

+ 2(< x X y 2 X xz > + <y X x 2 Xyz > + < z X xy > 2 + < x > 2 < y 2 z> + <y > 2 < x 2 z >) 
+ 4(< x X xy X yz > + < y X xy X xz > + < x X z X xy 2 > + < y X z X x 2 y >) 

I' 


+ 8 < x X y X xyz > +24 < x 2 y 2 z > 


S 

V , 


(10d) 


Expressions for integrals of other monomials in Cartesian coordinates may be obtained by appropriate 
substitution into the above formulas, and simplifying where possible. 
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III. High- Order Quadrature Approximations for the Fundamental Shapes 


Quadrature approximation for integrals of an arbitrary function f over the fundamental domains are of 
the form 


/ 


dL 

/( r) dS 

dV 


Y^ W 9^ T <l) 


L <7 



(11a) 


where r 9 is a quadrature point and the coefficient w q is the corresponding weight. The approximation has 
a degree of precision d if it is exact for all polynomials of degree equal to or less than d. Here we consider 
only approximations which are symmetric with respect to the vertices. The quadrature points then fall into 
symmetry groups, each of which is associated with a single weight. The quadrature approximation then 
takes the form 

dL 

f(r)dS = 
dV 


/■ 



( 116 ) 


- 9 Q 

where r| is a quadrature point belonging to symmetry group g. The symmetry is most clearly evident 
when viewed from the centroid of each fundamental shape. We therefore find it useful to introduce a 
parameterization based on the centroid, rather than the conventional simplex coordinates, to classify the 
quadrature points. 

Let r denote the position vector with respect to the centroid of the fundamental shape. Therefore 


where 


r - r - r c , 

(12) 

c - 1 E r * 

X 

(13) 

M 

►ii 

II 

o 

(14) 


defines the centroid. It is easy to show that 


From the above equation, it follows that in evaluating the integrals of tensor products of r , all summations 
in which any index appears once also vanish. For the line segment, since fi — — r 2 , we then obtain 

J fdL ~ j r r r dL = J fffrr dL = 0. (15a) 

The remaining eqs. (6) and (7) in terms of centroid based position vectors reduce to 


/ 


f- 

dS 

r 

dV 

r 

dL 

/ rr 

dS 

J 

dV 

f— 

dS 

/ rrr 

dV 


dL 

f rrrr 

dS 

j' rrrr dV 


dS 

rrrrr 

dV 


= 0 



;I>Ws) 

i x ' 

- — [V y' (f iTifjfj + •••) + 6 V] f ,] v 

n(n+l)(n + 2)(n + 3) l ^Y ^ 

48 ^ ( S ^ 

n(n + l)(n + 2)(n + 3)(n + 4) 4- W' 


(156) 

(15c) 

(15d) 

(15e) 

( 15 /) 

( 155 ) 
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Note that in general the integrals of the tensor products of f up to the fifth order are just proportional to 
the sum over all vertices of the tensor products of r The only exception is for the integral of the fourth 
order tensor product over the tetrahedron, where the additional of terms 10 ,j(.*i*i*j*j + * * *) (total of 
three terms) are necessary. 

A. Symmetry Groups 

The first symmetry group consists of the centroid r c . The second symmetry group consists of the n 
one-parameter vertex based quadrature points 

if = af t ( “ <Q<1 ,q^0). (16) 

a = 1 corresponds to the vertex, while a = — ^ry ^ opposite vertex, midpoint of the opposite edge, and 
centroid of the opposite face for the line segment, triangle, and tetrahedron, respectively, a = 0 corresponds 
to the centroid, and therefore it must be excluded from this group. Note that for the line segment (n = 2), 
the parameter a and —a describe the same group. If we sum the contribution from all members of the vertex 
based group, we obtain, from eq. (1), the relation 

E/( f n=n/(r c ) + y[^r,f ! ]:(VV/) c + ^[y]f £ f,f i ]:(VVV/) c + (17) 

i i i 

The third symmetry group consists of the n(n - 1) two-parameter edge based quadrature points 

f if = 7ft +6fj ( j ^ i, 7 £ <5 ^ 0) (18) 

for edge ij. To exclude points which are already in the first two groups, 7 and 6 must not equal zero, and j 
must not equal i. Here we also assume that 7 and 6 are not equal. (When the two parameters are the same 
we defined a fourth group, to be discussed below.) For the line segment, the third group is not necessary, 
since the first two groups cover the complete domain. It is easy to show that the values of 7 and 6 in the 
7 — 6 plane are restricted to the triangle determined by the points (0, 1), (1,0), and (-— ^ » —^2) in order 
for the quadrature points to lie in the domain. When 7 4* 6 = 1 the quadrature points lie on the edge ij. 
Note that for the triangle (n = 3) the group (7, 6 ) is also equivalent to the groups (-(5, 7 — 6) and (6 - 7, -7). 
If we substitute eq. (18) into eq. (1) (expanded about the centroid) and sum over all the edges, we obtain a 
series involving sums of tensor products of fT 5 . These can be expressed in terms of sums of tensor products 
of ii as 


t jjti 


EE f 7 « 

l j^l 


V V f^r^r 7 * 

tj tj tj 

l j^l 


EE 


-7<5-7$-7<5-7<5 
ij ij ij ij 


i j*i 


XT' V 

ij ij ij ij ij 

i j^i 


K 7 2 + <5 2 ) ~ (7 + <§) 2 ] r ' r ' 

i 


[n( 7 3 + «5 3 ) - (7 + <5) 3 ] * iTiTi 

i 

' [n( 7 4 + <5 4 ) - (7 + <5) 4 + 127 2 6 2 ] £7 fjfifjf, (n=3) 
(n(7 4 + 6 4 ) - (7 + <5) 4 ] 

< Z (n=4) 

2 7 2 ^ 2 + ' ' •) 

^ i j 

[n(7 5 + <5 5 ) - (7 + ^) 5 + 12(7 3 <5 2 + 7 2 <5 3 )]^f t f t r i f i f i . 


(19a) 

(196) 

(19c) 


(19d) 

(19e) 
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In deriving the above equations, we made use of the same simplifications that were employed to obtain 
eqs. (7). The fourth symmetry group is a degenerate case of the third symmetry group when 7 = <5 = /?, 
and consists of the n(n — l)/2 one- parameter edge based quadrature points 

ffj = P(fi + fj) (j>i, -P- ( 20 ) 


Q = \ corresponds to the midpoint of the edge ij. This group is not necessary for the triangle (n = 3), since 
the vertex based group already includes this case. We therefore show results for the tetrahedron (n = 4) only. 
For this case the parameter (3 and —(3 describe the same group. The relevant tensor product summations 
over the edges can then be expressed as 


EE'S-* 

i j>i 

i j>i i 

i j>i 

E E •S'S'S'S - ME E< f < va + - ■> -‘E™ 


i j>i 

i j * ij * t j * i j * ij 


z j 


EE^.«.-»' 


i j>i 


(21a) 

( 216 ) 

(21c) 

(21d) 

(21c) 


The conventional simplex coordinate /a, are defined by 


where 


r = E ^ ri ’ 

i 

E* = L 

i 


(22g) 

(226). 


It is easy to establish the relations between our centroidal parameterization and the conventional simplex 
coordinates for the four symmetry groups as follows: 


For 


f 76 


Vi = - 
n 


(for all i) 


Vi = 


1 + (n - l)a 


Vj = 


1 - a 


(j 7 ^ 0 


Mi 


Vj 


1 + (n - 1)7 - 6 


n 


1 — 7 + (n — 1)<5 


Vk — 


n 

I—7-6 


n 


(fc / i, fc ± j) 


4 : 


Hi = Hi = 


1 + (n - 2)0 


Hk = 


1 - 2(3 


n 


(k^i,k^ j) 


(23) 

(24a) 

(246) 

(25a) 

(256) 

(25c) 

(26a) 

(266) 
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The results contained in eqs. (15), (17), (19), and (21) can be summarized in table 1. It shows coefficients 
multiplying the various tensor product summations for each of the symmetry groups, as well as the integral 
of /( r) over the fundamental domains. Note that for n = 2 or 3, the column + •**) is 

not applicable, but six times the coefficients in that column must be added to the coefficients in column 
r if if if i. All symmetric quadrature formulas of any degree of precision up to five can be derived from this 
table. The number of independent terms (columns) required for a given degree of precision depends on the 
fundamental shape. For the line element, the degree of precision d of the quadrature approximation equals 
one plus two times the number of independent terms. For the other two shapes, an approximation of degree 
of precision d has d independent terms, except for the tetrahedron, where we have d + 1 terms if d > 4. 

In a Gauss quadrature formula, the parameters defining a symmetry group are unspecified. Then the 
number of unknowns associated with a group is one (the weight) plus the number of parameters defining the 
group. The number of quadrature points belonging to a symmetry group is a function of the fundamental 
shape. The efficiency of a particular symmetry group in forming a Gauss quadrature approximation is 
measured by the number of quadrature points per unknown. These are tabulated for the four symmetry 
groups and three fundamental shapes in table 2. For the line element, the two symmetry groups are equally 
efficient, but for the other two shapes the efficiency is a function of the symmetry group, with the centroid 
being the most efficient and the two-parameter edged based group being the least efficient. The latter group 
is not required to derive the most efficient Gauss quadrature formulas when the degree of precision is no 
greater than five. But that group is important in obtaining efficient Lobatto formulas, as will be shown 
below. Using tables 1 and 2 we can derive Gauss quadrature formulas for the three fundamental shapes by 
equating the total number of unknowns to the number of independent terms. There may be more than one 
combination of symmetry groups satisfying this condition. The combination giving the most efficient formula 
will be referred to as the Gauss quadrature formula. Since the parameters occur nonlinearly for d > 2, there 
may be more than one solution for a particular combination of symmetry groups. This will be found to be 
true for one case even after discarding solutions that locate quadrature points outside the domains. For the 
triangle and tetrahedron, the formulas may contain negative weight coefficients. 

While the Gauss quadrature formulas are the most efficient for isolated fundamental shapes, they may 
not be appropriate if positive weight coefficients are desired or if the shapes are part of a grid in which their 
boundaries (vertices, edges, or faces) are shared by several neighbors. In these cases it may be more efficient to 
employ limiting forms of the symmetry groups in which the points lie on the boundary elements. Quadrature 
formulas using these special symmetry groups will be referred to as Lobatto. These new symmetry groups 
consist of the vertices r V, edge midpoints rf^, face centroids r{, and one-parameter edge points r^. They 
are defined as 


and 


= rf with a = 1; 

r^ = r£ with a = — vertex k opposite edge ij, for triangle; 
with /3 = for tetrahedron; 


■«6 


f s = f° 
1 % ~1 


with a = — i, vertex i opposite face /, for tetrahedron; 


(27) 

(28a) 

(286) 

(29) 


r 4 4 = r7? with 7 = g, 6 = 1 - e. 




1; 


(30) 


In order to determine the relative efficiency of these new symmetry groups, we need the relations among the 
fundamental elements (vertices, edges, faces, and cells) which make up an arbitrary polygon or polyhedron. 
For a polygonal face (or plane polygon in two dimensions), which can always be triangulated without 
introducing interior vertices, we define 


v = number of vertices. 


(31a) 



e = total number of edges 
e\y = number of boundary edges 


(316) 

(31c) 


/ = total number of triangular facets or faces 


These satisfy the following relations: 


v = ej> 
e — 2eb - 3 


/ = e b - 2. (32c) 

For polyhedra, we restrict ourselves to triangulated polyhedra that can be subdivided into tetrahedra without 
introducing interior vertices. This is true for convex polyhedra. The subdivision may require introducing 
interior edges. For this case, in addition to the elements defined in eqs. (31), we also define 

e* = number of interior edges (33a) 

fly = number of boundary triangular faces (336) 

and 

c = number of tetrahedra (33c) 

The elements defined by eqs. (31) and (33) satisfy the following relations for a triangulated polyhedron: 

v = J + 2 (34a) 

e = c + f b + 1 (346) 

e h = \h (34c) 

6i = e — efc (34d) 

and 

/ = 2c+^. (34c) 

J 2 

The efficiency of the new symmetry groups in a Lobatto quadrature approximation is also a function 
of the type of integral being approximated. There are three such integrals in a conservation law. If the 
basis functions are not readily integrable, a quadrature over the conservation domain is required in step 3 
of the general discretization procedure. Since the same basis functions are utilized for all the cells in the 
grids, we consider integration of arbitrary functions over all the tetrahedra (or triangles) contained in the 
tetrahedralized (or triangulated) grid. Neglecting the effects of grid boundaries, and assuming a homogeneous 
grid, we obtain for domain integrals in tw T o dimensions the relations 

v = h (35a) 


e ”f >■ 

In three dimensions, one obtains analogously the relations 

1 

v = -c 
6 

7 

e= 6 C 



/ = 2c. 


(36c) 


The second type of integral is a source integral. This involves integrating an arbitrary non-linear function 
over a polyhedral (or polygonal) cell. The relations among the geometric elements are given by eqs. (32) 
and (34). The third type of integral is a boundary integral. For each function evaluation at a boundary 
point, several additional projections onto the boundary normals or tangents may be required. We assume 
that the computational cost of these projections can be neglected compared to the computational cost of 
the function evaluation. For two-dimensional problems the boundaries are line segments, and the Gauss 
formula is to be preferred. For three-dimensional problems, the required relations for the boundary of a 
triangulated polyhedra are given by eqs. (34b) and (34c). The data for each boundary symmetry group 
obtained from eqs. (32), (34), (35), and (36) is tabulated for the three types of Lobatto integrals in table 3. 
For completeness, we have also repeated the information for the group r c and rf . All Lobatto formulas can 
be obtained using tables 1 and 3, where entries in table 1 are replaced by their limiting values as indicated 
in eqs. (27) to (30). For a given degree of precision, there are now many more combinations of symmetry 
groups to evaluate as possible candidates. Again, we may find several or no solutions, and will discard those 
with quadrature points outside the domains. 

The last column in table 3 deserves further comments. For a given polyhedron, increasing the number 
of interior edges increases the total number of tetrahedral cells, faces, and edges, but the total number of 
vertices remains the same. Therefore, the net result is to increase the total total number of quadrature 
points for every symmetry group except the group r?, Whenever possible, one should triangulate polygonal 
faces and subdivide the triangulated polyhedron in such a way as to minimize the number of interior edges, 
which also minimize the number of tetrahedral cells. We illustrate this point with some simple examples. A 
tetragonal pyramid with a planar base, which has five vertices, can only be subdivided into two tetrahedra 
with no interior edge. But a trigonal bipyramid, which also has five vertices, can be subdivided into two 
tetrahedra with no interior edge, or three tetrahedra with one interior edge. The former subdivision is to be 
preferred, since it results in fewer total number of quadrature points. There are two important shapes with 
six vertices. An octahedron can only be subdivided into four tetrahedra with one interior edge. A trigonal 
prism can be subdivided into three or four tetrahedra, depending on the triangulation of the quadrilateral 
faces. It is easy to prove that an arbitrary trigonal prismatic grid can always be triangulated so as to result 
in three tetrahedra per prism with no interior edge. Another example is a structured grid consisting of 
hexahedra with quadrilateral faces. As shown by Rizzi (ref. 7), there exists a triangulation of the faces 
which permits subdivision into five tetrahedra with no interior edge. One of the five is twice the size of the 
other four. Kordulla and Vinokur (ref. 8) recommended a subdivision into six tetrahedra, by introducing 
a main diagonal as a common interior edge. Although the latter produces a more uniform subdivision, the 
former is preferred. In summary, polyhedra consisting of less than four tetrahedra will not require an interior 
edge, while those with four or more may require one. 

B. Quadrature Formulas 

With the aid of tables 1-3, we can derive virtually all symmetric quadrature formulas of any degree 
of precision up to five. By choosing different combinations of symmetry groups one might end up with 
more unknowns than the number of equations required to be solved. The resulting system of equations thus 
becomes under-determined, and the solutions can then be expressed in terms of free parameters. We can 
show that many published formulas are just solutions with particular choices of the free parameters. We have 
tabulated some Gauss and Lobatto formulas that are useful for approximating the different types of integrals 
for arbitrary polyhedral and polygonal grids. Only those formulas for which the quadrature points lie within 
the domains are given. The formulas for the line element are classical, and therefore are not presented. The 
coefficients for triangles and tetrahedra are given in tables 4 and 5, respectively. Although they can all be 
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expressed in closed form, due to space limitations, one of the relations is given in the text below. The letter 
g refers to Gauss formulas, while the other letters indicate Lobatto formulas. Those formulas not found in 
references cited in this paper are starred. Note that negative weights are present in formulas 3 g and 4c for 
triangles, and 3<?, 4 g, 46, 4c, and 5a for tetrahedra. 

We now give brief derivations of some the formulas, concentrating on those combinations of symmetry 
groups resulting in the solution of quadratic or cubic algebraic equations. As can be seen from tables 2 and 
3, the edge-based symmetry groups r T 5 and rf^, and their limiting forms are less efficient than the centroid r c 
and vertex-based symmetry group rf , and its limiting form. None of the edge-based groups or their limiting 
forms are required until we reach d > 4. We note from table 1 that if the centroid r c is present, the weight 
w c appears only in the first column. The system is thus simplified, having been reduced by one equation 
and one unknown. A unique solution is then readily obtained for all the cases we tabulated except for d = 5, 
where two formulas require the solution of a quadratic equation. 


It follows from table 1 that for d < 3 the triangle and tetrahedron may be treated in a unified manner. 
From table 2 we see that an n- point Gauss formula for d = 2 is given by the symmetry group rf . FYom the 
first two columns of table 1 we obtain immediately 


d = 2: 


w a 


l 

n 


and a = ± 

y/n 4* 1 


(37) 


For the tetrahedron (n = 4) the negative root must be discarded, since a is outside the permitted range (see 
eq. (16)). For the triangle (n = 3) both roots are valid, but the negative root (a = —\) places the points on 
the mid-edge of the triangle. According to our definition it is still a Gauss formula, with a Lobatto flavor, 
since it actually employs the mid-edge symmetry group r4. It therefore has superior properties, and is to 
be preferred to the positive root. From table 3 we see that it is also superior to the true Lobatto formula 
based on r c and for the source and boundary integrals, but equally as efficient for the domain integral. 
Using tables 1 and 2, we readily derive the solution for the n - 1-1 point Gauss formula for d = 3, based on 
r c and rf, as 


3 : 


W c = - 


n(n 4 l) 1 


w a = 


(n + 2 ) 2 
4 n(n 4- 1) 1 


Q = 


n + 2 


(38) 


Note that the centroid has a negative weight for both shapes. The Lobatto formula based on and rf 
involves the solution of a quadratic equation. The final expressions are 


1 =F >/4n 4 9 2n 2 4 6 n 4 3 ± >/4n 4 9 _ 1 ± \/4n 4 9 . 

^ Wy 2n(n 4 l)(n 4 2) ’ Woc 2n(n 4 l)(n 4 2) 7 2(n 4 2) 

The two roots for a are within the permitted range for both shapes, but only the lower sign gives all positive 
weights. For the tetrahedron, there is an additional reason for choosing the lower sign, since we have a = — ^ , 
so that the formula is actually based on rf and r{, and all the points lie on the boundary. On the other 
hand, one can show from table 3 that for the triangle the Lobatto formula based on r c , rf, and rfj is more 
efficient than eq. (39) for the domain and boundary integrals. 


When d > 4, it follows from column 5 in table 1 that the formula for a tetrahedron requires the symmetry 
groups rfj or fT 5 , or their limiting forms. The use of the two-parameter edge-based group is in general less 
efficient than using the combination of one-parameter edge-based and vertex-based groups. When employing 
F^-, one can derive from table 1 the relation 



( 40 ) 
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One immediately obtains w e = ^ for Lobatto formulas involving the edge midpoint symmetry group f^. 
For these cases, the solutions for tetrahedra have the same structure as those for the corresponding triangles. 

All the quadrature formulas except two for d > 4 that we have tabulated include two vertex- based 
symmetry group r^ 1 and rf 2 , or their limiting forms. Here we present a general method for obtaining 
solutions, although many Lobatto formulas can be determined by simpler procedures. The equations to be 
solved are some or all of the set 


w ai + w a2 = ao 

(41a) 

w ai Qi 2 + w a 2 q 2 2 = a 2 

(416) 

w ai Qi 3 + w Q2 a 2 3 =a 3 

(41c) 

WaiCtl 4 + l«a 2 a 2 4 = a 4 

(41d) 

u; Ql ai 5 + u; a2 Q 2 5 = a 5 . 

(41®) 


The parameters a* could include unknowns from other symmetry groups. From eqs. (41b) and (41c) we 
obtain for w ai and w a2 the expressions 


a 2 Q2 — O' 3 a 2 G:i — &3 

^ct i * o/ \* » ^<*2 o / T ■ 

ai 2 (a 2 -ai) a 2 2 (a i-a 2 ) 

Substituting the above equations into eqs. (41a), (41d), and (41e) we obtain 

a 3 (Qi + a 2 ) - a 2 [(ai + a 2 ) 2 - aia 2 )] = a 0 (aia 2 ) 2 
a 3 (ai "1“ <22) — a 2 aia 2 = 04 
a 3 [(ai + a 2 ) 2 - aia 2 ] - a 2 aia 2 (ai -j- a 2 ) = a 5 . 


(42) 


(43a) 

(436) 

(43c) 


For d = 4, we can eliminate (c*i + a 2 ) by substituting eq. (43b) into eq. (43a), and obtain for aia 2 the 
quadratic equation 


{a 2 z - aoa 3 2 )(axa 2 ) 2 + 2a 2 (a 2 a 4 - a 3 2 )a ia 2 + a 4 (a 2 a 4 - a 3 2 ) = 0. (44) 

If there axe no additional unknowns from other symmetry groups, and for the values of ao to a 4 in our 
tabulated cases, eq. (44) has two real solutions. For each of them, eq. (43a) then yields a quadratic equation 
whose solutions are ai and a 2 . In our cases the constants are such that one of the solutions of eq. (44) 
results in locating one a outside the permitted range. For d = 5, it follows from eqs. (43b) and (43c) that 
a\ and a 2 are two solutions of the quadratic equation 


(a 2 a 4 - a 3 2 )a 2 + (a 3 a 4 - a 2 a 5 )a + (a 3 a 5 - a 4 2 ) = 0. 


(45) 


As an example, we present the Gaussian solution for the tetrahedron for the degree of precision five. The 
most efficient quadrature points for this case involve two vertex-based r f and one edge-based symmetry 
groups. (The Gauss formula given in ref. 5 is less efficient than the present one.) We must now include 
eq. (40) in our equation set (41). In terms of 



(46) 


the parameters a t become 


_ 1 _ _A2_ 1 A 

a ° — 4 560’ fl2 - 20 420’ 

1 1 ^ 

03 = 77;, a i = Z 7 > and fl 5 = 


1 

140' 


(47) 
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Substituting (ai + 02 ) and a\Oc 2 from eq. (45) into eq. (43a), one obtains a cubic equation for A 


9A 3 - 284A 2 + 2800A - 8512 = 0. 


( 48 ) 


The above equation has three real roots. The values of a\ and c *2 from one root are complex, while those 
from another root place one set of points outside the tetrahedron. The remaining root is given 


A = 


27 


4\/79 


cos 


f cos- 1 67v ^ +2 x N 
COb 24964 ^ Z7r 


+ 71 


(49) 


This value must be used in evaluating formula 5# in table 5. Note that this formula was derived in reference 2 
by an entirely different procedure, but only numerical values of the parameters were given. We have provided 
explicit exact relations for the parameters. 


For a given degree of precision d, the optimum quadrature formulas to use (in terms of the minimum 
number of quadrature points) depend on the type of integral being approximated, and whether negative 
weights are allowed. For each value of d, using tables 1-3, we have examined all possible candidate combi- 
nations of symmetry groups to determine the optimum one. The results are summarized in table 6. In the 
w column, “pos” means that positive weights are necessary, while “neg” means that negative weights are 
allowed. For triangles, the formulas for the domain integrals were previously given in reference 4. One final 
note concerns the domain integral for a tetrahedron. For a non-uniform grid obtained by subdividing the 
hexahedron of a structured grid into five tetrahedra, the entries in table 3 for the groups r}\ r^, and 
would be |, |, and -y, respectively. This does not affect the choice of optimum formulas in table 6. 


IV. Concluding Remarks 


We have presented formulas for exact integrals of polynomials and quadrature approximations up to 
degree five for the line segment, triangle, and tetrahedron. These three shapes are the building blocks for 
integrations over arbitrary polyhedral or polygonal grids. A table of coefficients for quadrature symmetry 
groups and integrals for all three shapes is given, enabling one to construct symmetric quadrature formulas 
in a rational manner. The new procedure for obtaining the solutions in higher dimensions is no more complex 
than the one we are familiar with for one dimension. We have tabulated many useful Gauss and Lobatto 
formulas for both the triangle and tetrahedron, including a number that have not appeared in references 
available to us. All the formulas derived in this paper provide necessary tools in the spatial discretization of 
finite- volume equations over arbitrary polyhedral or polygonal grids to an accuracy of up to the fifth order. 
Higher-order formulas can be readily obtained using our procedure. Some of the quadrature formulas will 
now require the introduction of the two-parameter edge based symmetry group defined by eq. (18). As a 
result, the coefficients will not always be expressible in closed form, but may require the numerical solution 
of non-linear algebraic equations. 
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(c) tetrahedron 
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Figure 1. Parametric representations for the three fundamental shapes. 
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Table 1. Coefficients for quadrature symmetry groups and integrals for the fundamental shapes. 
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Table 3. Number of unknowns and Lobatto quadrature points for symmetry groups. 
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Table 4. Quadrature formulas for a triangle. 
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Table 5. Quadrature formulas for a tetrahedron. 
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Table 6. Optimum quadrature formulas for various types of integrals. 
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